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Abs tract 

The so-called H-type mesh is used in a finite-element (or finite- volume) 
calculation of the potential flow past an airfoil. Due to coordinate 
singularity at the leading edge, a special singular trial function is used for 
the elements nei^boring the leading edge. The results using the special 
singular elements are compared to those using the regular elements. It is 
found that the unreasonable pressure distribution obtained by the latter is 
removed by the embedding of the singular element. Suggestions to extend the 
present method to transonic cases are given. 
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1. Introduction 

In finite-volume calculations of the potential flow around an airfoil 
(Jameson and Caughey, 1977), the physical space is mapped to a computational 
space in which the airfoil is mapped to a horizontal or vertical line segment 
and a rectangular mesh is created. The velocity potential in the rectangular 
element is then represented by a bilinear function of the computational 
variable for a two-dimensional calculation. This process assumes that the 
transformation from physical space to computational space is regular every- 
where. Because of this assunption, the type of mesh one can use is restricted. 
For an airfoil calculation, a C-type or 0-type mesh which wraps around the 
leading edge is commonly used, while ^e H-type mesh (see Figure 1) which has 
a transformation singularity at the leading edge is usually considered not 
suitable for calculation. hkDwever, the latter has some advantages in the 
geometrical description of a complex three-dimensional configuration. An 
obvious example is the canard-wing cootoination. In order to place the shed 
vortex on the prescribed surface, the wakes coming off the canard and the wing 
must be on mesh lines. It is impossitile for a C-type or an 0-type mesh to 
fulfill this requirement. For an H-type mesh, this requirement can be satis- 
fied fairly easily. Since we know that the difficulty with an H-type mesh 
stems from the singularity of the geometrical transformation, there must be a 
way to "undo" this singularity. 

In this report, we represent the potential and the transformation function 
near the transformation singularity by a trial function which reflects the 
local singular behavior in the computational space. This singular behavior 
cannot be represented by the conventional bilinear trial function and has been 
the heart of the problem in using the H-type mesh. Since the problem we are 
dealing with is common in all Mach number ranges, we have chosen to work with 
inconpressible flow. The method developed here can, however, be extended to 
the transonic case. 



- 2 - 


Flow Research Report No. 216 
January 1982 



Figure 1. H-type Mesh 
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2. Finite-Element Formulation 

Consider the two-dimensional Laplace equation 

-<t> - d> = 0 (1) 

^xx ^yy 


with the far field Dirichlet boundary condition at infinity and the Neumann 
boundary condition on the airfoil. Variational principle for this problem 
states that solving Equation (1) together with the boundary conditions is 
equivalent to minimizing the quadratic functional 


Kf) 



dx dy 


( 2 ) 


for appropriate functions f which are square integrable [i.e., the right-hand 
side of Equation (2) is finite] together with their first-order derivatives 
(Strang and Fix, 1973). 

In the following finite-element formulation, we first generate a mesh of 
quadrilaterals in the physical space. For each cell w, the mesh is then mapped 
to a unit square fi in computational space by the transformation X = X(x, y) 
and Y = Y(x, y). It is more convenient to estimate 1(f) in ccmiputational 
space. Therefore, we shall need the following formula for the change of 
variable. Let f be any function on w and F be the corresponding function on 
fi, i.e., F(x, Y) = f(x, y). Applying the chain rule and the change of 
variable formula to Equation (2) we have 


Kf) 



(FyXx 



h"^ dX dY 


(3) 


where h is the Jacobian 


= Vv - 


( 4 ) 


The variational problem is now completely handled in the transformed 
space. The variation of F within the computational element Q is represented 
by some trial function having four corner values (o^, 012 , o^, a^) as para- 
meters. The elements are divided into two groups: regular elements in which 
the function F can be locally approximated by a bilinear function and special 
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elements In which F as well as the mapping (x, y) -»> (X, Y) are singular and 
special functional forms must be chosen. 

In the following, we shall deal with the interior elements first. 
Treatment of the boundary elements is given in Section 3. 

Regular Elements 

A schematic representation of the transformation between the physical 
space and the computational space is shown in the following figure. 


PHYSICAL COMPUTATIONAL 




Figure 2. Transformation Between Physical 
and Computational Space 


The local relationship between physical and computational coordinates is 
X = I - X) ^ - y) x^ + ^ + X) I - y) X2 

+ - X) + Y) X 3 + ^ + X) + Y) x^ . (5) 

A similar expression holds for y. 

In general, a bilinear function on the square cell (1, 2, 3, 4) with corner 
values , 02 » the form 

F = ^ - X) - Y) o^ + + X) - Y) 02 

+ (j-X) (“ + Y)o3 + ^+X) ^ + Y) 


“4 


( 6 ) 
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Using Equation (5), the expression for h in Equation (4) can be written as 
a linear function in X and Y: 

h = ho + hj^X + hyV (7) 

where , 

^0=7 ^^2 “ ^1 ’'a “ ^3^ (^3 ■ yi ^4 " ^2^ 

“ "2 (^3 ” ^1 ^4 ” ^^2 ~ ^1 * ^4 ” ^ 3 ^ 

and 

hx = («2 - «i) (V4 - yj) - (*4 - Xj) (>2 - ( 5 , 

*\ = (^4 - ^2^ (V3 - Yi) - (xj - (Y4 - . 

The finite- volume scheme of Jameson and Caughey (1977) approximates h by 
its value at the center of the cell. Here, we try to extend the accuracy one 
order higher by approximating h“^ in Equation (3) with a linear function. 

This can be done easily by binomial expansion. 

Using Equations (5), (6) and (7), the integral of Equation (3) in a regular 
element can be obtained as a quadratic functional of the comer potential values 
0.21 Oj and a^. By summing the integrals from all the elements, 1(f) is approx- 
imated by a quadratic expression of the potentials at the mesh points (f’ij)* 
Minimizing I is equivalent to solving for a system of linear equations: 

If = 0 . (10) 

The coefficients of this set of linear equations form the so-called stiff- 
ness matrix. The general expression of this stiffness matrix is quite cumber- 
some and will not be given here. For an identity transformation, the weight 
in the stiffness matrix is distributed as shown in the following figure. 

It is the sum of a five-point formula ^d a rotated (45®) five-point formula. 


. 


-1 

-1 

-1 

-1 
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Figure 3. Stencil of Stiffness Matrix 
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Special Elements 

The transformation from physical to computational space is set up to be 
analytic in the four cells at the leading edge except at the leading edge 
itself, and it transforms these cells to four squares of unit area. 



It is well known that a function which is regular in the physical space 
behaves as in the transformed space. Thus, rather than using a bilinear 
function and packing the mesh densely at the leading edge to approximate the 
singular solution, we use four special elements. 

On the triangle (1, 2, 4) shown in Figure 5, consider the function 

F = + (c^2 ” + (a^ - a2)YX”^'^^ (11) 



Clearly such a function attains the values a^, and at corners 1, 2 
and 4; furthermore, it is linear on the edge (2, 4) v^hich would assure 
continuity between special and regular elements. Along edge (1, 4), two such 
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functions from neighboring triangles are identical so that continuity between 
special elements is preserved. Heuristically, should pick up r^'^^ 
behavior and the variation of corner values ci 2 ahd should reflect the 
angular variations. A similar formula with the reverse roles of X and Y 
can be applied to the triangle (1, 3, 4). 

Let f be the function in physical space corresponding to F, which is of the 
form of Equation (11). We now evaluate 1(f). By analyticity, we have the 
Cauchy-Riemann equations: 

^ “ Yy » ^ ” “Vx (12) 


Substituting these equations into the Jacobian expression given as 
Equation (4), Equation (3) can be written as 


1(f) 



(13) 


Note that Equation (13) is similar to Equation (2) except that we have F in the 
integral on the right-hand side instead of f. 

The four special elements are broken into eight triangles, in each of which 
the trial function is assumed to be of the form of Equation (8). A calcula- 
tion of the 3I/3a terms using Equatiai (13) for I gives the following 
expression for the contribution to the stiffness matrix from the element shown 
in Figure 5: 


. 

3ttj^ * 

“l 

1 

fn|<t 

1 

3 1 

4 03 + 2 «4 

31 . 

9 ^ 

38 

29 

^2 ‘ 

■12 “l 

^ I ?“2 

■17 “4 

li . 

30j • 

9 

-T 2 “1 

+ 

38 29 

I2 “3 “ I2 “4 

31 . 

1 

29 „ 

29 ^ ^ 52 „ 

'“4 ‘ 

2 “l 

■ I 2 “2 ‘ 

I2 “3 ^ T2 “4 


We can take advantage of the permutational symmetry of the special elements to 
obtain the contribution to the stiffness matrix by the other three. 

The element stiffness matrix for each cell is computed and the results are 
assembled to give the global stiffness matrix K. 
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3. Boundary Conditions 

Let 4 be a potential which satisfies Equation (1) together with the boundary 
conditions. Let = x cos a + y sin o be the mdisturbed potential with flow 
angle a. Instead of solving for 4, it is more convenient to solve for the 
reduced potential G = $ - 4^. The boundary conditions on 4 must be transformed 
to those on the reduced potential G. 

Far Field Boundary Condition 

The far field boundary condition is given Iteratively as the potential 
vortex field generated by the lift. Consider the far field boundary cell 
(1, 2, 3, 4) shown in Figure 6. 


BOUNDARY 


3 


1 


4 


2 


Figure 6. Far Field Boundary Cell 

Let edge (3, 4) be part of the boundary. In solving for the solution 
values, which are comer values of G, and are given by the potential 
vortex obtained from the previous iteration. Hence, one assembles the stiffness 
matrix for this boundary cell in the same way as for an interior cell. In com* 
puting the solution, terms involving and are known quantities and are 
moved to the right-hand side of the equation. 

Neumann Boundary Condition 

The Neumann boundary condition on the airfoil for 4 is 4 = 0, where n 
is the outward normal unit vector to the airfoil (See Figure 7). While the 
boundary condition for 4 is the so-called natural boundary condition in the 
variational method, that on G requires some modification. The corresponding 
boundary condition for G is 


3^ G = - 3;^ 4„ 
n n ® 


( 15 ) 




Flow Research Report No. 216 
January 1982 


-9- 



Figure 7. Cell With Neumann Boundary Condition 


Consider now a cell with the Neunann boundary condition on edge (1, 2). 

Then, Equation (15) together with the expression for takes the form 

3-j- G = - n, cos a - n« sin a (16) 

n 1 2 

where n^^ and r\^ are x- and y-components of n [the outward normal of the line 
segment (1, 2)]. 

For a cell with the Neunann boundary condition as in Figure 7, we define 

Kg) = f * ^y)^^ ^ W + '^2 ^ 

cell *' 


(17) 


where ^da is the integration on the line segnent (1, 2). We next reproduce 
the well-known proof that if G is a function which minimizes Kg), then G 
satisfies Equation (16) (Strang and Fix, 1973). Let c be a real number, then 

KG + eH) = + Gyj dx dy + dx dy 

+ 2e j + GyHy) dx dy + 2-^[n^ cos a + Oj sin ct) (G + cH) ctr 

= 1(G) * 2eJ + GyHy) dx dy + 2c cos ot + n 2 sin a)H da 

+ + H^Jdx dy . (18) 

Since G minimizes Kg), the e term in Equation (18) must vanish, i.e., 

/ [=x'^x * «yV '<>' * COS a +02 sin a)H da =0 (19) 
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for any function H such that 1(H) < *. Integrating by parts, we have 


+ -j- G) + (nj^ cos a + nj sin a) 


H da = 0. 


( 20 ) 


Since H is artitrary except for the constraint that 1(H) < «, each integral 
in the above expression must vanish independently. The vanishing of the 
second Implies Equation (16). 

Thus, cells with Neumann boundary conditions differ from interior cells 
only by the second term on the right-hand side of Equation (17). Using the 
linear expression for G on the line segment (1, 2), this term can be estimated 


as follows: 


2 (n^^ cos a * ^2 sin a) G da 

rl/2 

= 2L / (nj^ cos a + n 2 sin a) - X) + (^ + X) 02 

4/2 


dx 


= L (n^ cos a + n 2 sin a) + 02 ) 


( 21 ) 


where the first equality comes from Equation (6) , and the second equality 
comes from a straightforward calculation. L is the length of the line segment 
( 1 , 2 ). 

For the two special elements with Neumann boundary conditions, we obtain 
results similar to Equation (21) by using a highest order approximation. On 
the boundary of the special element cell, the analytic mapping can be approxi- 
mated by C = as shown in Figure 8. 



Figure 8. Local Mapping of Boundary Near Singularity 
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Then, 


1 - 1/2 

da = dC = LX dX 


Using this together with Equation (11), we have 
2 ^('^1 ^2 “) G dd 


= L (n^ cos a + 02 sin a) I (“2 ” °1^ X^^^ 

O*' ^ 

= L cos a + 02 sin a) + 02 ) 


X‘^^^ dX 


The result is identical to Equation (21). 

Thus, for cells with Neumann boundary conditions, the calculation of 
Sl/Sa^ gives an extra term of the form 

L (oj^ cos a + 02 sin a] . 


( 22 ) 


(23) 


( 24 ) 


When solving for the values, this term is known and is moved to the 
right-hand side of the equations. 
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A. Mesh Generation 

The mesh is generated by the following procedure. The curvature ic of the 
airfoil at the leading edge is estimated by using the leading edge and two 
adjacent data points. A parabola of the same curvature through the leading 
edge is fitted in the physical space z = x 4 iy. Next, the parabola is mapped 
onto a semi-infinite line (0, s) in C space (C = s + it) by the inverse of 
the analytic function 

z = i + IC C . (25) 

We can now generate a mesh of rectangles in the % space such that the four 
cells at the leading edge form four squares. The mesh in C space is mapped 
back to z space by relation (23). Finally, we use a shearing transformation 
to shear the parabola to the airfoil surface to obtain a body-fitted mesh 
system (see Figure 9). 



Figure 9. H-type Medium Mesh 
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5. Numerical Calculations and Conclusions 

After gathering the stlffhess matrix, the solution to the llr^ar equation 
system Is obtained by a horizontal line relaxation scheme. The relaxation 
starts at the line extending from the leading edge to the upstream Infinity. 
Since the leading edge Is a mesh singularity, the mesh line extends downstream 
from the leading edge along both the upper side and the lower side of the 
branch cut, l.e., along the upper and lower surfaces of the airfoil 
(Figure 10). The potential values along this entire line are updated simul- 
taneously by a special solver. Then, the ordinary trldlagonal solver Is used 
to relax the lines successively away from the airfoil. 



In order to verify the accuracy of the code, a test case for flow past a 
NACA 0012 airfoil Is carried out with an angle of attack of 3 degrees and 300 
Iterations for each of the coarse, medium and fine meshes. The results (see 
Figures 11-16) match those of FLO 26 for a C-type mesh (finite-volume scheme, 
Jameson and Caughey, 1977). For the calculation using an H-type mesh with 
leading edge cells treated as regular elements, the pressure distribution 
obtained exhibits an unreasonable "kink" near the leading edge as seen In 
Figures 12 and 13. This "kink" seems to persist as the mesh Is refined. 
Figures 14 and 15 show the corresponding results using special elements In 
four cells neighboring the leading edge. These results show that the un- 
reasonable behavior of the pressure distribution near the leading edge Is 
removed. For further comparison we also Include the result by a panel method 
(Figure 16). Again, the results agree extremely well with those obtained by 
the finite-element method with singularity embedding. 
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Figure 11, Finite-Volume Scheme (FLO 26) 
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Figure 12. Regular Eletnent Method, Medium Mesh 
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CL 0.3636 CD 0.0051 CM 


-0. 0046 


Figure 13. Regular Element Method, Fine Mesh 
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Figure 14. Special Element Method, Medium Mesh 
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Figure 15. Special Element Method, Fine Mesh 
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Figure 16. Panel Method 
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A few words can be given to the possibility of extending the present method 
to transonic cases. The best way seems to be adapting the artificial density 
concept. The conservation of mass can be written as 

V • (p V<j>) = 0 (26) 

where "p = p - yp^^ (27) 

is the artificial density and pgAS is the upwind difference of the 
isentropic density. The nonlinear equation is solved iteratively as 

V • =0 (28) 

where the superscript n denotes the nth approximate solution. For the (n+l)th 
approximation, is a known function. The variational principle is then to 
minimize 

= / (p V(|> • V(|)) dx dy . (29) 

0) 

To the degree of accuracy comparable to the existing finite- volume scheme, 

"p in an element can be approximated by its value at the center of the cell. 

The singularity embedding method can be extended to the transonic case by 
multiplying the incompressible cell stiffness matrix by its cell artificial 
density. One also needs to modify Equation (17) for Neumann boundary condi- 
tions. The rest of the method remains unchanged. One can try to obtain a 
more accurate result by approximating the isentropic density by a linear 
function in a cell and treating the added part (-yp^AS) as the difference in 
the center value. The stiffness matrix will be more complicated, but no 
essential difficulty is expected. 



Flow Research Report No. 216 
January 1982 


- 21 - 

References 

Jameson, A., and Caughey, D. A. (1977) "A Finite Volume Method for Transonic 
Potential Flow Calculations,” in Proceedings of the AIAA 3rd Computational 
Fluid Dynamics Conference. Albuquerque, June, pp. 35-54 (AIAA Paper No. 
77-635). 

Strang, G., and Fix, G. J. (1973) An Analysis of the Finite Element Method. 
Prentice-Hall . 


WJ:589R 



1 Report No. 

NASA CR 166387 

2 Govamment Accession No 

3 Reapient's Catalog No 

1 4 Title and Subtitle 

Singularity Embedding Method in 
Potential Flow Calculations 

1 

5 Report Date 

June 1982 

6 Performing Organization Code 

7 Author(s) 

Wen-Huei Jou and Hung Huynh 

8 Performing Organization Report No 

10 Work Unit No 

T-9682 

11 Contract or Grant No 

A-77471B 

13 Type of Report and Period Covered 

Contractor Reoort 

14 Sponsoring Agency Code 

9 Performing Organization Name and Addrtta 

Flow Research Company 
21414-68th Ave. South 
Kent, Washington 98031 

12 Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 


]S, $uppi(tntnt*n Notes 

Point of Contact: Dochan Kwak, NASA Ames Research Center, Moffett Field, CA 

94035, (415) 965-6743, FTS 448-6743 


16 Abstract 

The so-called H-type mesh is used in a finite-element (or finite- 
volume) calculation of the potential flow past an airfoil. Due to 
coordinate singularity at the leading edge, a special singular trial func- 
tion is used for the elements neighboring the leading edge. The results 
using the special singular elements are compared to those using the regular 
elements. It is found that the unreasonable pressure distribution obtained 
by the latter is removed by the embedding of the singular element. Sugges- 
tions to extend the present method to transonic cases are given. 


17 Key Words (Suggested by Author(sl) 

Singularity Method 
H-grid 

Computational Aerodynamics 

18 Distribution Statement 

Unclassified - Unlimited 
STAR Category 02 

19 Security Qassif (of this report) 

Unclassified 

20 Security CUssif (of this page) | 

Unclassified 1 

21 No of Pages 

24 

22 Price* 

. 


*For sale by the National Technical information Service. Springfield Virginia 22161 






















End of Document 



